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Abstract 

We propose a minimal model of the dynamics of diversity — replicator equations 
with extinction, invasion and mutation. We numerically study the behavior of this 
simple model and show that it displays completely different behavior from the con- 
ventional replicator equation and the generalized Lotka-Volterra equation. We reach 
f-^ , several significant conclusions as follows: (1) a complex ecosystem can emerge when 

mutants with respect to species-specific interaction are introduced; (2) such an 
ecosystem possesses strong resistance to invasion; (3) a typical fixation process of 
{SJ . mutants is realized through the rapid growth of a group of mutualistic mutants 

with higher fitness than majority species; (4) a hierarchical taxonomic structure 
(like family-genus-species) emerges; and (5) the relative abundance of species ex- 
hibits a typical pattern widely observed in nature. Several implications of these 
results are discussed in connection with the relationship of the present model to the 
generalized Lotka-Volterra equation. 

X" 

5— i 

Key words: replicator equation, generalized Lotka-Volterra equation, diversity 



1 Introduction 



The relationship between the complexity and stability of an ecosystem has 
been one of the most fascinating topics in theoretical biology for decades 
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(Pimm 1991). In the 1950s and 1960s the proposition that highly complex com- 
munities are more stable than simple ones was widely supported (MacArther, 
1955; Elton, 1958). However, this view was challenged by theorists in the 1970s, 
who discussed the stability of a community of species interacting randomly 
(Gardner & Ashby, 1970; May, 1972, 1974; Roberts, 1974; Gilpin & Case, 
1976; Taylor, 1988a). The most significant result was given by May (1972), 
who considered a large-dimensional ecological equation with an n-dimensional 
random interaction matrix {%} whose diagonal elements are —1 and whose 
off-diagonal elements are assigned as Gaussian random numbers with mean 
and variance a 2 (with probability C) or (with probability 1 — C) . He found a 
kind of phase transition from stability to instability in which the equilibrium 
solution corresponding to the coexistence of all species becomes unstable if 
a > {nC)~ l l 2 . He concluded that an ecological system cannot be stable if it 
is complex. After May's work and the clear conclusions he reached, a number 
of mathematical biologists attempted to explain the discrepancy between the 
observed complexity of ecosystems in nature and the results of these mathe- 
matical studies. Their works can be classified into two groups with regarded 
to the methodology of assembling a complex and stable ecosystem. 

The first pioneering idea was formulated by Tregonning & Roberts (1979). 
They prepared a 50 x 50 random interspecies interaction matrix and checked 
the feasibility of the equilibrium solution of a population equation. Here, an 
equilibrium solution is termed feasible if all species have positive populations. 
In the case that the solution under consideration was not feasible, they re- 
moved that species whose equilibrium population was most negative. Repeat- 
edly applying this procedure, they acquired a stable equilibrium of positive 
population with considerably higher diversity than that predicted by May's 
theorem [see also Roberts & Tregonning (1980)]. Let us call their approach 
the eliminating approach. Although this approach can produce a complex and 
stable ecosystem, the resulting diversity is in fact much less than the initial 
diversity (= 50). Moreover, the biological significance of their modeling was 
not very clear because negative population is never achieved in nature and be- 
cause the origin of the random interactions was not discussed. [See also Tokita 
& Yasutomi (1999) for a paleontological explanation of the formation of such 
random interspecies interactions.] Furthermore, according to the analysis of 
extinction dynamics (Tokita & Yasutomi, 1999), which is a type of population 
dynamics similar to the elimination process, the resulting diversity is indepen- 
dent of the initial diversity. This result suggests that when the interactions 
are assigned randomly, a system with arbitrarily large diversity cannot emerge 
through such an elimination process. However it is still worthy of attention 
that Tregonning & Roberts (1979) reported that the most often-observed inter- 
species interaction of the survived subsystem was a prey-predator relationship, 
because it was recently demonstrated that coexistence of a number of species 
can be achieved by an ecological model with such prey-predator relationships 
even if the interactions are random (Chawanya & Tokita, 2002). 
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Other authors have used a second method, which is known as community as- 
sembly (Pimm 1991). In this method, a model community is assembled through 
a sequence of invasions of species. This type of invasion is believed to reflect 
the natural tendency of species to arrive sequentially rather than simulta- 
neously. The method further can be classified into two subgroups. The first 
trailblazing one is represented by that of Post & Pimm (1983). They prepare 
a pool of species ex ante and then add them sequentially. The properties of 
the newly introduced species are then determined through a sampling process 
that is completely independent of the properties of the species existing in the 
community at the time of introduction. A considerable number of studies have 
followed this line (Robinson & Valentine, 1979; Yodzis, 1988; Taylor, 1988b; 
Drake, 1990; Case, 1990; Law & Morton, 1996; Happel & Stadler, 1998). These 
models are realistic in the case that we consider systems in which new species 
come from outside of the existing community. We thus refer to this as the 
invasion approach. 

The second subgroup concerns a longer-term evolutionary process (Colwell & 
Winkler, 1984; Ginzburg et al, 1988; Happel & Stadler, 1998). In such models, 
a new species is introduced as a mutant of an existing species. More precisely, 
parameters characterizing a new species are generated by making appropri- 
ate changes to the parameters of one of the existing species, rather than by 
sampling the parameter values independently from a separate predetermined 
ensemble. We refer to this method as the mutation approach. 

Although these works of community assembly succeed in producing stable 
communities whose numbers of species are larger than that predicted by May's 
theorem, the coexistence of hundreds of species has not been demonstrated. 
On the other hand, it is worthy of mention that in a recent mutation approach 
(Drossel et ai, 2001), considerably higher diversity has been realized than that 
by other mutation approaches. They focused on influence of prey-predator 
interactions to the evolution of food web structure. Consequently, mutualisms 
and other interspecies relationships were not considered, and "connectance" 
(the percentage of non-zero values of interactions) was low. In contrast, in 
the present study, we consider not only the prey-predator type but also high 
connectance of interactions of mutualisms, competitions and parasitisms. 

In the present model, species are defined by interactions matrices, following 
the preceding studies mentioned above. Such definition is, however, still a 
controversial paradigm because interaction coefficients are not real biological 
traits, like body size, but rather summaries of these traits filtered through 
some interaction rules. We neverthless adopt the interaction matrices because 
one of the motivations here is to clarify the relationship between complexity 
and stability of large ecosystems, i.e. the paradox of ecology, raised by May 
(1972). For that purpose, analyses of multidimensional trait space would be 
essential, which is, in general, not a simple task for the real biological traits. 
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Here we only cite some approaches on community assembly that does deal in 
real traits (Brown & Vincent, 1992; Sasaki 1997; Geritz, et al, 1998). 

In this paper, we observe the behavior of a system based on the replicator 
equation in which a species is eliminated when its population becomes too 
small, and mutants or invaders are added periodically. We consider the cases 
of three different rules governing the generation of new species. We call these 
invasive, global and local. The invasive rule can be compared with the model of 
the invasion approach, and the global rule with the global mutation approach. 
The local rule is introduced here for the first time. 

We will show that the local rule exclusively allows for the assembly of a very 
large and complex ecosystem which contains hundreds of species strongly in- 
teracting with each other. In addition, we observe that the diversity and mu- 
tualism increase together. The key point in our model is the introduction of a 
mutant species whose relationships with the dominant species (with large pop- 
ulations) are virtually the same as those of their parent species, but who can, 
nevertheless, interact differently to minor species (with small populations) as 
well as to those species which will come in the future. The fitness of such a 
mutant is almost identical to that of its parent when it is introduced, and 
hence it often behaves as a temporally neutral mutant. 

We will show that the minor mutants play a key role in the emergence of 
a complex ecosystem. This idea can be discussed in connection with neutral 
molecular evolution (Kimura, 1983), which stresses the importance of the ac- 
cumulation of nondirectional mutations. Our theory is closely related to the 
ecological neutral evolution theory for the abundance and diversity of species 
in tropical rain forests and coral reefs proposed by Hubbell (1979, 1997, 2001). 
The size distribution of the relative abundance of species obtained in our model 
is in good agreement with a well-known distribution obtained from detailed 
observations of real ecosystems. 

We will also discuss the resistance to invasion of assembled ecosystems. We 
find that an ecosystem assembled under the local rule has perfect resistance 
to invasion if we keep it free from invaders which are assembled according 
to the invasive rules for a sufficiently long time. Contrastingly, an ecosystem 
assembled under the global rule cannot resist invaders. 

This paper is organized as follows. In the next section we observe the dy- 
namics of the number of surviving species in our model. In the third section, 
we discuss the resistance to invasion of the assembled ecosystems. The final 
section is devoted to discussion of the relationship between the generalized 
Lotka-Volterra equations and the replicator equations. 
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2 The replicator equation with extinction and mutation 



2.1 Model 



Here, we investigate a model based on the system of ordinary differential 
equations 



dxi(t) 



AT/ Nj 



Xi(t) \ Y. a ii x A t ) - 2 ajkXj(t)x k (t) ] , 



called the replicator equations (RE) (Hofbauer & Sigmund, 1998), on a JV/- 
dimensional simplex: 

AT/ 

= 1 (o<^(o<i). (2) 

i=l 



The variable denotes the population (or relative abundance) of the species 
i, and Afj denotes the initial number of species, that is, the initial value of the 
diversity. The (i, j)-th element of the matrix A = (a^) determines the effect 
of species j on the growth rate of species i. 

It is known that the A/-dimensional RE is topologically equivalent to an 
Nj — 1 dimensional generalized Lotka-Volterra equation (GLVE) (Hofbauer & 
Sigmund, 1998). We use the RE here, simply because they are more tractable 
for a numerical simulation of dynamics of a large ecosystem with high diversity 
and complex interactions. A comparison between the GLVE and RE in the 
context of the dynamics of diversity is given in the last section. 

It should be noted that the RE may possess heteroclinic orbit even in low 
dimensions {Nj > 4) (May & Leonard 1975; Chawanya 1995 & 1997). When 
a heteroclinic orbit approaches a saddle, where some species are 'extinct', 
their population take extremely small values but they never exactly reach 
zero, because the orbit is bound in the interior of the simplex (2). In the 
vicinity of such a saddle, these population have such small values that they 
cause underflow in naive numerical calculations. However, if we continue a 
calculation in such a region using the more sophisticated calculation method 
of Chawanya (1995, 1997), we often observe some of these species start to 
revive and cause the orbit to eventually leave for another saddle, in particular 
for a system with high diversity and complex interactions. This transition 
among saddles continues cyclically or chaotically. The exponential approach 
of a population to zero and its subsequent revival to order 0(1) plays an 
important role in heteroclinic orbits. However, in the real world, such small 
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population cannot be realized. 



Considering this problem, we introduce a parameter 5 1), the extinction 
threshold, into the dynamics described by (1) and (2) : At each discrete time 
step, the population Xk is set to zero if this quantity becomes less than 5. 
The population of the surviving species {x{\ (i ^ k) are then renormalized to 
satisfy J2i=ik x i = 1- The diversity N(t) decreases through the above described 
process. We refer to this as the extinction dynamics (ED) model (Tokita & 
Yasutomi, 1999). With 5, a stochastic effect on extinctions is introduced into 
the fully deterministic RE, and hence, 5 plays a role analogous to that of 
random genetic drift. The introduction of S constitutes a finite size effect on the 
total population into RE, because S represents a minimum unit of reproduction 
for each species, and its reciprocal 1/8 represents to the permissible population 
size. That is, the system explicitly possesses an energy constraint and is free 
from the problem of population explosion often observed in the simulation of 
a large-scale GLVE pointed out by Taylor (1988b). 

In contrast to the typical ED simulation with a large Nj, in the present sim- 
ulation, we start from only one species (i.e. Nj = 1), following Ginzburg 
et a/.(1988). The initial intraspecies interactions is set to — v (< 0), which is 
one of the parameters in the present model. A simulation proceeds by first 
computing the ED for a fixed number (T) of time steps. This is followed by 
the introduction of a new species. This procedure is then repeated many times. 
We set T = 900 in our simulations. Hereafter, we refer to T time steps as one 
period, and we use r to count the number of periods. At the end of each period, 
we select one parent species, j, among the existing species [1 < j < N(r)] with 
a sampling probability proportional to the population Xj. Then we produce a 
new species k [= N(r) + 1)] from the parent species j, copying and {cir,} 
f° { a ki} and {an-} for all i, with some variations applied in the manner de- 
scribed below. We then set the new species' population as x k = ( (( > 8) and 
normalize all Xi in order to maintain the condition J2 x i — 1- Thus completing 
the introduction of the new species, we then return to the ED, and the entire 
process is repeated. 

As we discussed above, we consider three different rules for the variations of 
the interactions {aki} and {a^} for the new species: invasive, global and lo- 
cal. In the case of the invasive rule, we add a new species whose interaction 
coefficients are assigned randomly and are independent of all existing species 
(including its parent). In the case of the global rule, we produce a mutant whose 
interaction coefficients are all only slightly different from those of its parent 
species. In the case of the local rule, we introduce a mutant that is identical 
to its parent species, apart from its interaction with one particular surviving 
species. While the invasive (Robinson & Valentine, 1979; Yodzis, 1988; Happel 
& Stadler, 1998) and global (Ginzburg et al, 1988; Happel & Stadler, 1998; 
Drossel et al, 2001) rules have been employed in previous studies, the local 
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Fig. 1. Only the local mutation rule allows the number of observable species to 
continue to increase, (a) The horizontal axis r represents the time measured with 
respect to the interval between 'mutations'. The vertical axis represents the number 
of observable species S(t). (b) The log-log plot of S(t) under the local rule from (a) 
and the function with which it is fitted, r ' 38 (dashed curve). The parameters used 
here are 5 = exp(— 11), £ = exp(— 10), £ = exp(— 9), v = 1.0 and T = 900. 



rule is examined for the first time in the present study. In terms of theoreti- 
cal population genetics, the invasive rule corresponds to the "house-of-cards" 
approximation (Turelli, 1983). The details of these rules are as follows. 

• Invasive: We assign the values of {aki} and {a^} for all i (7^ k) as Gaussian 
random numbers (mean and variance 1). This implies that the mutant 
species is independent of its parent species. 

• global: Quantities {vki} and {vik} representing Gaussian random noise with 
mean and variance 1 are independently added to {a^} and {a^} for all i 
according to 

' a ki <^= ^a ki + (1 - ^)v ki 
\aik <= la-ik + (1 - l)vik 

for all i, where 7 (0 < 7 < 1) denotes the strength of the noise. In our 
simulations, this value was set to 0.9. 

• local: We select a species / whose population is non-zero and change only 
the relationship between the mutant species and species / (7^ k). Thus, the 
effects of a mutation are concentrated on only two elements, a k i and a^. We 
assigned the values of these two elements using Gaussian random numbers 
with mean and variance 1. 

In all cases, the intraspecies interaction an takes the value — v for each species 
i. Finally, we introduce a threshold £ (£ > ( > 5), and we consider a species 
to be unobservable if its population Xi is smaller than £. Since £ is larger than 
C, a mutant species is not observable until its population increases and be- 
comes greater than £. The function S(t) denotes the number of the observable 
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Fig. 2. (a)The average fitness, f(t) = J2j t=i a jkXj(t)xk(t). All curves were obtained 
from an average over 50 samples of independent simulations with different series of 
random numbers, (b) The log-normal plot of / under the local rule from (a) and 
the function with which it is fitted, log(r ' 25 ) (dashed line). 

species existing at the end of the rth period. We regard S(t) as the diversity 
in this model. It should be noted that, in general, there appear much more 
unobservable species than observable species although the population of the 
latter is much larger than that of the former. 



2.2 Results 



Figure 1 (a) displays the development of S{t) for each of the three rules, 
where S{t) is defined as S{r) averaged over independent simulations with 
different series of random numbers. We set the initial conditions as X\ — 1 
[that is, iV(0) = 5*(0) = 1] and a n = —v. This figure shows that diversity 
grows most under local rule. The diversity S{r) under the local rule increases 
to approximately 50 at r = 10, 000, and as can be inferred from this plot, 
it continues to increase beyond this value as more mutants are added. We 
averaged the diversity over 50 simulations and found that this averaged value 
obeys S(r) ~ r 38 [Figure 1(b)], which should be compared with the square 
root law S(t) ~ y/r in the coinfection model (May & Nowak, 1995). 

Figure 2 (a) displays the development of the average fitness under the three 
rules. As in the case of the diversity, it is seen here that only the local rule 
allows for the average fitness to continue to increase. Since the average fitness 
/ = Y,fk=i ajk x j( T )%k(T) is a kind of average of the interaction coefficients, we 
can regard this quantity as an index of mutualism. The results here therefore 
indicate that the diversity and the level of the mutualism increase together. 
It is interesting that the average fitness does not display power-law behavior, 
but rather increases much more slowly: / ~ clog(r) (c ~ 0.25) (Figure 2(b)). 
This is analogous to the slow increase of the diversity observed in a GLVE 
model of super- infection in a host-parasite system (May & Nowak, 1994). 
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Fig. 3. (a) Development of the diversity S(t) and the average fitness /. This figure 
shows the result for a single sample under the local rule, (b) An enlarged view of 
(a) for r = 7500 - 7900. Here 5 = exp(-ll), ( = exp(-lO), f = exp(-9), v = 1.0 
and T = 900. 

Figure 3(a) plots the development of S(t), where the largest number of species 
observed for this case is 196 at r = 144, 700. In the plots of Figures 3 (a) 
and (b), only one sample for the development of each S(t) and / is used. 
By contrast, the plots in Figure 1 represent averages over 50 samples. From 
Figure 3(a), we wee that the diversity does not increase monotonically, and 
repeated avalanches of extinction followed by rapid recoveries are observed. 
Figure 3(b) gives an enlarged view of the time evolution plotted in Figure 
3(a) from r = 7,500 to 7,900 in order to clarify the relationship between 
S(t) and /. Here we observe a clear negative correlation; that is, the average 
fitness drastically increases during an avalanche of extinction. This behavior 
is typically observed in extinction dynamics (Tokita & Yasutomi, 1999). How- 
ever, the average fitness does not increase during the time that the diversity 
increases. This implies that extinction plays a more important role than spe- 
ciation in the increase of the average fitness. This is because large extinction 
events swiftly exterminate less mutualistic species and produce a more mu- 
tualistic network of species. The refined network then prepares for a rapid 
recovery of the diversity in the aftermath of the extinction. 

In the case that S(t) is denned with just one sample, its fluctuation suggests 
simultaneous occurence of extinction and speciation of multiple species. On 
the other hand, /(r) develops relatively smoothly, because the definition of 
/ is a kind of average of the interactions {<Zy} over {xi}, and Xi at the time 
that the ith species appears or extincts is very small (8 for extinction and ( 
for speciation). 

Figure 4 plots the dependence of S on the ratio of the value of the diagonal 
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Fig. 4. The observable diversity S(t) at 
r = 10, 000 displayed as a function of 
v under the local rule, v represents the 
value of diagonal elements of the in- 
teraction matrix (the intraspecies in- 
teraction), while the off-diagonal in- 
terspecies interactions are assigned as 
Gaussian random numbers with mean 
and variance 1. Here S = exp(— 11), 
£ = exp(— 10), £ = exp(— 9) and 
T = 100. Each dot corresponds to an 
average over 10 samples of indepen- 
dent simulations. 



elements, v, and the variance of the Gaussian random numbers (= 1) that 
are used to assign the interaction coefficients of the mutant in the case of the 
local rule. We see that S depends only weakly on v. We should note that May 
(1972) showed that the internal equilibrium point becomes unstable when this 
value becomes less than 1 if the system is assembled randomly at one time, 
rather than through the periodic introduction of mutants, as done in our 
study. The strength of the interactions between species affects the stability of 
the ecosystem if it is assembled at one time, but it does not do so when the 
ecosystem gradually develops for a long time on the evolutionary time scale. 

Why does the diversity S(t) increase so drastically only in the case of the 
local rule? To answer this question, let us consider Figure 5, which gives an 
example of our numerical integration. This figure shows that a few species 
enter the existing ecosystem more rapidly than exponentially and virtually 
simultaneously. This kind of rapid population growth results from the fact 
that the relationship between these mutant species is more mutualistic than 
those between the dominant species. When the relationship between mutants 
p and q is more mutualistic than those between the dominant species, and this 
relationship therefore makes these mutants "fitter" than the average species, 
their population increase. The growth of the population of species p thus in- 
creases the fitness of species q, thereby enhancing its growth rate, and hence 
increasing the fitness of p. These unobservable species constitute a new mutual- 
istic unobservable subsystem and eventually merge into the observable existing 
ecosystem. The implication of this behavior is that the formation of a mutual- 
istic network is more effective as a invasion strategy than competition or even 
exploitation under the pressure of natural selection, because a mutualistic in- 
vasion has the abovementioned positive feedback effect, while the exploitation 
involves a negative feedback (as the more one species preys on another, the 
less abundant this prey becomes). 

In order to form such a subsystem, mutants must wait for the arrival of part- 
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Fig. 5. A few mutant species grow 
more rapidly than exponentially, and 
they enter into the observable ecosys- 
tem at virtually the same time. Here 
5 = exp(-ll), C = exp(-lO), 
f = exp(-9), v = 1.0 and T = 900. 
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Fig. 6. Mutant species produced through 
the mutational rule grow less rapidly 
than exponentially. After entering the 
ecosystem, they usually come to replace 
their parent species. Here S = exp(— 11), 
C = exp(-10), £ = exp(-9), v = 1.0 and 
T = 900. 



ners. In order to avoid extinction while waiting, such a mutant must have 
virtually the same fitness as its parent species. More precisely, they must have 
the same interaction coefficients with respect to the dominant species as their 
parents. If this is not the case, they will independently invade the existing 
ecosystem, replacing their parent species, or they will simply be repulsed. Of 
course, the invasive and global rules rarely produce such "sleeping" mutants, 
because in general in these cases all interactions of a given mutant will be 
different from those of their parents. Figure 6 gives an example of our nu- 
merical integration for the case of the global rule. Here a new species invades 
independently, replacing an existing species, and we can observe no faster- 
than-exponential population growth. Only the local rule can produce such 
mutants, and this is the reason why only it can yield a remarkable increase in 
the diversity. 

Since a subsystem of the type described above has a higher degree of mu- 
tualism than the observable system, we expect that the ecosystem will be- 
come more mutualistic as the diversity increases. Since the average fitness 
1 aj] c Xj(t)xk(t) is a kind of average of the interaction coefficients, 
we can regard this average as an index of mutualism. As seen in Figure 1, 
it continues to increase only under the local rule. Figure 7 exhibits the dis- 
tribution of {%•} for an ecosystem that has emerged under the local rule. It 
should be noted that these values are 'naturally selected' from Gaussian ran- 
dom numbers with mean and variance 1. As seen in the figure, the average 
of the distribution shifts to the positive area. This suggests that the relation- 
ships between species tend to become more and more mutualistic through the 
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Fig. 7. The vertical axis represents the 
frequency, and the horizontal axis repre- 
sents the values of the off-diagonal ele- 
ments of the matrix (a^) produced un- 
der the local rule. The distribution of 
the elements of this interaction matrix 
shifts toward the positive direction as 
time increases. The distribution is av- 
eraged over 10 samples of simulations. 
Here 5 = exp(— 11), ( = exp(— 10), 
£ = exp(-9), v = 1.0 and T = 900. 



2 - 



1 - 




• Local 




Gaussian 

I ' 1 1 1 1 

12 3 

V", 

Fig. 8. The distributions of the for 
each species under invasive, global and 
local rules at r = 10, 000. For reference, 
the (described in the main text) 
also appear. As seen here, the under 
the local rule deviate greatly from the 
fW. Here 5 = exp(-ll), £ = exp(— 10), 
£ = exp(-9), v = 1.0 and T = 900. 



struggle for survival. 

Let us consider in detail the process of speciation by introducing a 2S(r)- 
dimensional trait space, where each species i is characterized by the point 
(aii, Oi2, • • • , a»,s(r)) a Ui a 2i, ■ ■ ■ , a s(r),i)- In order to visualize the distribution of 
the species in this high- dimensional space, we project this point onto a point 
in a two-dimensional space: 

^>4<>)^fe ^4) (4) 



If the ecological status of two species k and I are similar, b^ k ' and are close 
to each other. 

Figure 8 displays the sets {b^} for each species with the interaction matrices 
resulting from invasive, global and local rules. For the sake of contrast, we 
also display the points denned in the same way as the &W, for a 50 x 50 
matrix whose off-diagonal elements are given as Gaussian random numbers 
with mean and variance 1. We can see that the sets {b^} for these three 
interaction matrices differ greatly from {fW}. This suggests that it is quite 
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Fig. 9. (a): The genealogy of surviving species at r = 50,000 and 100,000. Each 
species is plotted as a point whose location is determined by the time that this 
species emerged (horizontal axis) and the generation that this species represents 
(vertical line). Mother species and their daughter species are connected by lines. 
The connections between existing species at r = 50, 000 and their ancestors appear 
as black lines, and those at r = 100,000 as gray lines, (b): This figure shows the 
following; (1) species closely related in lineage possess similar characteristics (at 
least to the extent that these are captured by the quantity plotted here); (2) there 
is a hierarchical structure reminiscent of a family-genus-species; and (3) almost all 
species existing at r = 50, 000 are extinct at r = 100, 000, and the descendants of 
only three species survive and produce a large number of descendents (adaptive radi- 
ation). Black circles indicate {frW} at r = 50, 000 and white circles at r = 100, 000. 
Here S = exp(-ll), ( = exp(-10), f = exp(-9), v = 1.0 and T = 900. Species 
which belong to a given branch are enclosed within the same closed curve. The thin 
closed curves correspond to new branches, and the thick closed curves correspond 
to old branches. Thus the thin curves enclose a single genus and the thick curves 
enclose a family. 

inappropriate to use a random matrix to approximate an interaction matrix of 
a real ecosystem. Among these three sets {b^}, that obtained from the local 
rule in particular is shifted in relation to {r^} more mutualistic region in 
the trait space. This is suggested also by Figure 7. Moreover, several groups of 
species are observed in this case, indicating that the speciation processes under 
the local rule do not result in a simple random drift or a random diffusion in 
the trait space; that is, the groups of species that appear in this case are fixed 
by the pressure of natural selection. 

Figure 9 displays the relationship between species in the case of the local rule. 
It is observed that this ecosystem starts from a single species, and the diversity 
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Fig. 10. Rank-size plot of 

the population densities X{ 

for 178 observable species at 
r = 150, 000. 

50 100 150 200 

Rank 

increases through speciation. Figure 9(a) displays the history of the speciation. 
The black lines indicate a genealogy of species existing at r = 50, 000, and 
the gray lines at r = 100, 000 from their ancestors at r = and r = 50, 000, 
respectively. Lines representing the extinct species before r = 50, 000 and 
t = 100, 000 are not plotted, although there are in fact more branches for the 
extinct species. At r = 50, 000 there exist 58 species, and they are separated 
into two large groups which branches at a very early stage. 

Figure 9(b) displays the points for each species. The set {b^} at r = 
50, 000 is plotted with black dots and those at r = 100, 000 as circles. We 
observe two large groups in this figure: Groups A and B contain 39 and 19 
species at r = 50, 000, respectively. Those species on the same branch of 
the genealogy are plotted near each other in the 2-dimensional space. We also 
observe several subgroups in each group: Group A contains four subgroups and 
group B contains two. These groups also reflect the structure of genealogy, as 
those species in the same subgroup appear on the same sub-branch of the 
genealogy. If we think of each dot and circle as a 'species', then each subgroup 
can be thought of as a 'genus' and an entire group is a 'family'. The two 
families A and B still exist at r = 100, 000, shifting toward the upper right. 
All species at r = 50, 000 are extinct, and those descendants of only three 
species survive at r = 100, 000. Each arrow denotes a radiation from one of 
these three to a group or a subgroup at r = 100, 000. Descendants of one 
species of family A radiated into 32 species. This family separated into two 
genera, and even these genera formed subsubgroups. There are 13 species of 
descendants that emerged from two species of family B: one of these genera 
contains 4 species and the other 9 species. The latter genus is separated into 
4 sub-genera. 

Finally, we note the interesting correspondence of a result obtained in the 
present study to a well-known ecological law of population abundance: As 
seen from Figure 10, the dominance-diversity distribution is S'-shaped, and its 
middle part displays a power law behavior. Such behavior is observed in many 
species-rich communities (Hubbell, 1997; Magurran, 1988). This suggests that 
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the present model shares some characteristics with populations in nature. Let 
us here note that our data were not obtained through any averaging, but 
from only one simulation. This is significant because in order to obtain such 
a distribution using one sample, it is required that the sample produce a 
sufficiently large number of species. This has been realized for the first time 
in the present study. 

In the present study, we have mentioned "invasion type" mutants with respect 
to species-specific interaction by the local rule. However we can also think of 
"mutant type" mutants with respect to species-specific interaction, which have 
the parameter deviated from the parent's value by a random quantity. We have 
executed similar simulations by such a mutative local rule and found a similar 
diversification like the local rule. Under the mutative local rule, the diversity 
S(t) increases little slower but is fitted by a power function like under the 
local rule. 

From the above discussion, we conclude that an ecological system which is 
composed of many hierarchically structured species strongly interacting with 
each other can emerge through an evolutionary speciation process by the local 
rule. 



3 Resistance to Invasion 



3.1 Method 



As mentioned in the first section, Elton (1958) suggested that complex com- 
munities are more stable than simple ones. In this context, as pointed out by 
Case (1990, 1991), Elton's definition of the 'stability' is not the asymptotic 
stability of the assembled matrix but the resistance to invasion of communi- 
ties. We thus can understand his hypothesis as an assertion that species-rich 
ecosystems are more resistant to invasion by exotics than are species-poor 
ones. 

Case (1990, 1991) constructed models to check Elton's hypothesis from this 
point of view. He randomly assembled ecosystems and selected stable and 
feasible ones. Then he tested their resistance to invasion by adding new species 
generated by a rule very similar to our invasive rule. He concluded that a 
more diverse and more strongly connected ecosystem has greater resistance to 
invasion. 
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Fig. 11. After R periods of as- 
semblage under local and global 
rules, invaders generated by the 
invasive rule attack for 200 pe- 
riods. When R is small, the 
invaders succeed in destroying 
the existing ecosystem. However, 
When R is larger than 100, the 
ecosystem developed under lo- 
cal rule can maintain its diver- 
sity. Contrastingly, that devel- 
oped under the global rule ex- 
periences a reduction in diver- 
sity resulting from attacks by 
invaders. Here 5 = exp(— 11), 
C = exp(-lO), i = exp(-9), 
v = 1.0 and T = 900. 

3.2 Results 

We carried out a similar test on ecosystems formed under the global and local 
rules and compared the resistance to invasion of these ecosystems. There is an 
important difference between these two rules and the invasive rule. Under the 
global and local rules, new species are mutants of the existing species. However, 
under the invasive rule these are independent of the existing species, invading 
from without. In our tests, we added mutants generated by the local and 
global rules for only the first R periods, and then we let invaders generated 
by the invasive rule attack the ecosystem for 200 periods. R denotes a kind 
of maturity index of the ecosystem. 

The results of our experiment are depicted in Figure 11, which shows the 
dependence of the resulting diversity S(t) on R. As seen there, when R is 
small, the diversity is drastically reduced under either rule. However, when R 
is larger than 200, ecosystems developed under the local rule maintain their 
diversity, resisting invasion. Contrastingly, those developed under the global 
rule exhibit a reduction in their diversity for any value of R, as a result of 
destructive invasions. This result suggests that an ecosystem developed under 
the local rule has strong resistance to invasion if it is sufficiently mature. 

The difference between the local and global rules results from the difference 
in their levels of mutualisms. As we have seen in the previous section, an 
ecosystem evolving under the local rule is more mutualistic than one evolving 
under the global rule. We should note that the average fitness / here again 
turns out to play an important role as an index indicating the resistance to 
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invasion . If an ecosystem has an average fitness /, the fitness // of an invader 
must be larger than / in order for the invader to succeed. Since // can be 
approximated by a Gaussian random number with mean and variance 1, for 
an ecosystem with sufficiently large diversity, the probability for a successful 
invasion is given by the error function, 



Prob.(/ 7 > /) = Erfc(/) = 




which decreases exponentially as / increases. From Figure 1(b), it is possible 
to estimate the probability (5) by using the value of / for R ~ 200 at which 
the difference between the local and global becomes larger in Figure 11. The 
average fitness for the local rule, /(r = 200) ~ 1, gives Erfc(l) ~ 0.16, while 
that for the global rule, /(r = 200) ~ 0.25, gives Erfc(0.25) ~ 0.4. Figure 
1(b) indicates that / is a monotonically increasing function of R for the local 
rule, and hence we expect that a random invasion becomes almost impossible 
against a fully- matured ecosystem (i.e., one of large R) under the local rule. 
That is, only temporally neutral mutants can influence or create extinction 
events in an existing ecosystem evolved under the local rule. 



4 Discussion 



Generalized Lotka- Volterra equations vs. replicator equations 



The generalized Lotka- Volterra equations (GLVE) are more commonly used 
than the replicator equations (RE) in the field of mathematical biology. How- 
ever, it is known that an iV-dimensional RE is mathematically identical to the 
(N - l)-dimensional GLVE (Hofbauer & Sigmund, 1998) 

-f = vM + £ KVj) (i = 1, 2, • • • , N - 1), (6) 



where the ith species' population yi is defined in terms of the population in 
the RE (1), Xi by 

y^ = — ■ (7) 

X N 



The value of yi can take any non-negative real number. The interactions {bij} 
and the intrinsic growth (or decay) rate r\ here are defined in terms of inter- 
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Fig. 12. The vertical axis repre- 
sents the frequency, and the hor- 
izontal axis represents the values 
of the antisymmetricity Das(B) at 
r = 4,000,170,000 and 700,000 
and Das (A) at r = 4,000 
and 170,000 ("RE"). The value 
of the antisymmetry Das (A) a t 
r = 700, 000 becomes nearly zero. 
Here S = exp(— 11), ( = exp(— 10), 
£ = exp(— 9), v = 1.0 and 
5 T = 1,000. 



actions in the RE by 

bij = ciij — a,Nj ; (8) 
Ti = a iN - a NN =a m + v. (9) 

Any orbit of the GLVE can therefore be mapped onto an equivalent orbit of 
the RE, and any behavior exhibited by these two systems can be studied using 
either, at least in principle. 

Here we explain why we nevertheless study the RE rather than the GLVE. 
As described in the first section, Taylor (1988b) pointed out that most of the 
dynamics resulted in explosions of populations as proportion of species with 
Ti > increased or proportion of positive b^ increased. Moreover, the popula- 
tions of species with low fitness often take extremely small values that cause 
underflow in naive numerical computational scheme. Consequently, numerical 
analysis of the GLVE often encounters both divergence and underflow, in par- 
ticular for a system with high diversity and complex interactions. On the other 
hand, in the RE, population explosions are avoided by definition, because any 
orbit of the population {xj} is bound in the simplex J2i=i x i — 1- Moreover, 
in the present model, the problem of underflow is also avoided because the 
heteroclinic orbits, the cause of underflow, are excluded by the introduction 
of the extinction threshold (Tokita & Yasutomi 1999). These are the reasons 
we use the RE rather than the GLVE. 

Let us note that the meaning of the interaction coefficients a^- in the RE is 
different from that of the in the GLVE: even if the and take only 
positive values (mutualistic relationships), the by and bji can take negative 
values (competitive or exploitative relationships) depending on the values of 
ciij and a Nj in Eq. (8). As the growth rate of the species i is defined by the 
fitness J2j a ij x j deducted by its average J2j,k UjkXjXk in the definition of the 
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replicator equations (1), the ecological significance of the developed ecosystem 
in the present study should be discussed by the matrix B = (b^) in the 
GLVE than A = (a^) in the RE. We therefore study the nature of matrices 
B transformed from A by Eq. (8). Here we should note that there are N ways 
of transformation from A to B because of the arbitrary option of the species 
N in Eq. (7). To characterize a L x L matrix M = {m^}, here we define 
"antisymmetricity" as 

Das(M) s «^ -"*>% (10) 



The bracket ((■ • -))m = l(l-i) ^i=i ^>j&(' ' ') denotes the average over all off- 
diagonal elements of M. For the RE with interaction A at r, L becomes S(t), 
and L = S(t) — 1 for B of the corresponding GLVE. In general, the antisym- 
metricity Das{M) becomes 0, 1/2 and 1 when M is symmetric (m^- = rriji), 
randomly asymmetric (m^ ^ rriji) and antisymmetric (m^- = —rriji), respec- 
tively. Figure 12 plots the distribution of D AS {B) at r = 4, 000, 170, 000 and 700, 000, 
where 5(4,000) = 45,5(170,000) = 125 and 5(700,000) = 216. The value of 
the antisymmetry Das(A) of the RE at r = 4, 000, 170, 000 is also indicated. 
From the figure, the interaction matrix A of the RE turns out to approach 
the symmetric point (Das — 0; mutualism). On the other hand, the matrix 
B tends to approach the intermediate region (Das ~ 0.3) between symmetric 
and asymmetric. This implies that the system evolved under the local rule has 
complex ecological interspecies interactions including mutualism, competition, 
predation and parasitism. 

From the figure 7 and 12, we conclude that the frequency of the interspecies 
interactions of the emerged system can be approximated by a symmetric Gauss 
distribution with mean m(> 0) and variance 1 



1 / ( a .._ m )2\ 
p ( a ij) = -^ ex Py tJ ~2 J' ( for and a ij = a ji ) (11) 

a ii = -v(<0) (12) 

in the infinite limit of the number of species (N — > oo), where m > is 
the level of mutualism. This type of the replicator equation, in general, has 
a number of saturated fixed points, which is proved by the symmetry of the 
interaction (Hofbauer & Sigmund, 1998). 

Finally, let us see how {rj} and {b^} in the corresponding GLVE mutate in 
the local rule. First, values r^, {bn~} and {bki} for all % of a new species k are 
copied from rj, {bij} and {bji} of a parent species j, respectively. Second, by 
the definition of the local rule and the transformation (8) and (9), r k or b [k (and 
bki) for arbitrary I is replaced by a random quantity. Note that r k is replaced 
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by a random number with mean v and variance 1 while bit {bki) with mean 
and variance 1. From the transformation (8) and the distribution (11) and (12) 
of dij, the distribution of 6^ becomes a Gauss distribution with mean and, 
therefore, is not mutualistic. Minor mutualist species in the RE, therefore, are 
not necessarily mutualists in the GLVE, and do not necessarily lead population 
explosion. On the other hand, in Eq. (9) is biased by v(> 0) and the mean 
value of Ti increases over time because mean value of increases as indicated 
by Fig. 7. Accordingly, in terms of the GLVE, proportion of producers increases 
as time goes on in the local rule. As pointed out by Taylor (1988b), this may 
trigger population explosion. Another chance for explosion is implied in the 
transformation (7) when the iV-th "base" species goes extinct (xjv(£) — > 0), 
which breaks down the simulation of the corresponding GLVE. We stress here 
that even if the extinction of one "base" species causes the explosion in the 
GLVE, the evolution in the RE proceeds successfully and there are still other 
equivalent GLVE systems by transformations by other N — 1 "base" species as 
Hi = Xi/xM (M = 1,2, ■ ■ ■ , N — 1). The RE, therefore, traces a possible path 
of evolution avoiding the breakdown of the simulation which may occur in the 
GLVE. This methodological advantage of the RE is related to the circumstance 
that the dimension (the degree of freedom) of the RE is one higher than the 
GLVE and the total population is conserved (Eq. (2)). This suggests that 
such a conservative quantity, e.g. resource limitations, would be essential in 
the modeling of multispecies evolution using the GLVE. 
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